Milan population

In the following plot it is represented the evolution of the Milan population, according to the ISTAT census. The data can be found here. As we can see, the population between 1980 and 1990 decreased. In order to take into account the decrement of the number of people subjected to the risk, I computed the population day by day with a linear interpolation of the census dates and used this piece of data as an exposure in the models. With this information I computed the mortality rate as: \[ Mortality Rate = \frac{tot.mort}{population} \] Note that this is just a generic mortality rate, not a specific one. So, we are not considering the age distribution of the population. Given that within the period 1980-1990 the Milan population aged, even If the mortality at each age unchanged, the generic mortality rate would have increased.

Mortality

In the following plots we can see the trend of the mortality rate and the respiratory mortality rate through time. As we can see there is a strong seasonal effect, with a higher mortality in winter and a lower mortality in summer. We can spot some outliers in summer 1983 and in the first quarter 1986. In those periods the mortality was much higher than the one registered in the same period in other years.

In the following plot we can compare the mortality pattern through the year between different years. As we can see, they are quite similar.

Explantory variables

In the following plots we can see the trend of the explanatory variables through time. As we can see, all these variables have a strong dependency with time. That means that an eventual strong correlation between these variables and the mortality could be due to a spurious correlation. For example, during winter:

• people work more, so they are more stressed and the accidents on the job are more frequent;
• the light time is shorter, and that has an impact on the health;
• there is more traffic, so the car accidents are more frequent.

To deal with this problem, I inserted the period of the year in the model with a spline computed on the day of the year.

Observing the mean.temp, trend we can see that in summer 1983, it was hotter than other summers. That could partially explain the higher mortality of that period.

Considering interaction

I also investigated whether is better to consider or not the interaction between day.of.year and mean.temp in the spline component of the model (i.e. considering s(mean.temp) + s(day.of.year) or s(mean.temp, day.of.year)).

fit_gam_bayes_int <- stan_gamm4(tot_mort_prob ~ s(day.of.year, mean.temp) + factor(year) + factor(weekend),
                                family = gaussian, data = mort,
                                chains = CHAINS, iter = ITER, seed = SEED, cores = CORES)

As we can see, plotting the joint effect of mean.temp and day.of.year for the two models, they are quite similar.

To choose the better among the two, I used the Leave One Out Information Criterion (LOO-IC). According to it, the model without the interaction performs better.

loo_best <- loo(fit_gam_bayes_best)
loo_int <- loo(fit_gam_bayes_int)

loo_best

## 
## Computed from 1000 by 3652 log-likelihood matrix
## 
##          Estimate   SE
## elpd_loo -10241.1 46.7
## p_loo        32.4  1.7
## looic     20482.2 93.4
## ------
## Monte Carlo SE of elpd_loo is 0.2.
## 
## Pareto k diagnostic values:
##                          Count Pct.    Min. n_eff
## (-Inf, 0.5]   (good)     3650  99.9%   370       
##  (0.5, 0.7]   (ok)          2   0.1%   313       
##    (0.7, 1]   (bad)         0   0.0%   <NA>      
##    (1, Inf)   (very bad)    0   0.0%   <NA>      
## 
## All Pareto k estimates are ok (k < 0.7).
## See help('pareto-k-diagnostic') for details.

loo_int

## 
## Computed from 1000 by 3652 log-likelihood matrix
## 
##          Estimate   SE
## elpd_loo -10266.5 49.2
## p_loo        39.1  1.3
## looic     20532.9 98.3
## ------
## Monte Carlo SE of elpd_loo is 0.2.
## 
## Pareto k diagnostic values:
##                          Count Pct.    Min. n_eff
## (-Inf, 0.5]   (good)     3650  99.9%   439       
##  (0.5, 0.7]   (ok)          2   0.1%   510       
##    (0.7, 1]   (bad)         0   0.0%   <NA>      
##    (1, Inf)   (very bad)    0   0.0%   <NA>      
## 
## All Pareto k estimates are ok (k < 0.7).
## See help('pareto-k-diagnostic') for details.

compare(loo_best,
        loo_int)

## elpd_diff        se 
##     -25.4      12.6

Model checking

The following plots are aimed to check if the model hypothesis are verified in the data.

From the overlay plot and the statistics computed on the whole dataset the model looks fine.

Considering the statistics computed for each year, something strange appears. Looking to the plots for the standard deviation, we can see that in the years 1983 and 1986 the standard deviation in the data is much higher than expected, while in other years, like 1982, 1987 and 1988 it is lower than expected. That means that the model doesn’t describe properly the variance within years and the hypothesis of homoscedasticity is not verified.

From the following plot we can see that in the years 1983 and 1986 there are some outliers with a particularly high mortality. In 1983 this phenomenon happened during summer, while in 1986 it happened during the first quarter.